Here, I compare using raster to projected shapefile to check for differences.
library(raster)
library(sf)
library(ncdf4)
library(rmapshaper)
library(tidyverse)
library(diptest)
library(moments)
library(viridis) #colors
library(data.table)
library(hydroTSM) #hypsometric curves
library(gridExtra)
library(maptools)Checking rgeos availability: TRUElibrary(rgdal)
library(rgeos)
library(SpaDES)Paths set to:
options(
rasterTmpDir = '/private/var/folders/0m/zw4hywk9177dsnzhv4w0b5pm0000gn/T/Rtmp9C6ACH/raster'
reproducible.cachePath = '/private/var/folders/0m/zw4hywk9177dsnzhv4w0b5pm0000gn/T/Rtmp9C6ACH/Require/cache'
spades.inputPath = '/tmp/RtmpF3R1kl/SpaDES/inputs'
spades.outputPath = '/tmp/RtmpF3R1kl/SpaDES/outputs'
spades.modulePath = '/tmp/RtmpF3R1kl/SpaDES/modules'
)
loading reproducible 1.2.4
loading quickPlot 0.1.6
loading SpaDES.core 1.0.5
loading SpaDES.tools 0.3.6
loading SpaDES.addins 0.1.2
Default paths for SpaDES directories set to:
cachePath:
inputPath: /tmp/RtmpF3R1kl/SpaDES/inputs
modulePath: /tmp/RtmpF3R1kl/SpaDES/modules
outputPath: /tmp/RtmpF3R1kl/SpaDES/outputs
These can be changed using 'setPaths()'. See '?setPaths'.library(rnaturalearth)
library(rnaturalearthdata)
etopo_shelf_df <- readRDS("~/Documents/grad school/Rutgers/Repositories/shelf_habitat_distribution/etopo_shelf_df.rds")
#bring in bathymetry data frame for shelf regions
#LMEs
LME_spdf <- readOGR("LME66/LMEs66.shp") #spatial points data frame with all 66 LMEsOGR data source with driver: ESRI Shapefile
Source: "/Users/zoekitchel/Documents/grad school/Rutgers/Repositories/shelf_habitat_distribution/LME66/LMEs66.shp", layer: "LMEs66"
with 66 features
It has 9 fields
Integer64 fields read as strings: OBJECTID #convert to equal area projection
#The Lambert azimuthal equal-area projection is a particular mapping from a sphere to a disk. It accurately represents area in all regions of the sphere, but it does not accurately represent angles.
equalareaprojection<- crs(" +proj=laea ")DF to Raster for bathymetry
etopo_shelf_raster <- rasterFromXYZ(etopo_shelf_df, crs = crs(LME_spdf))
#reclassify all values <2000m in depth to 1 instead of actual depth
etopo_shelf_raster<- reclassify(etopo_shelf_raster,cbind(-Inf, Inf, 1))
We will test regions that are likely to give us issues because they include high latitudes.
##Degree Shifts
Eastern Atlantic
LMEs to include -19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 58, 59, 60, 62
east_atl <- c(19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 58, 59, 60, 62)
east_atl_spdf <- LME_spdf[(LME_spdf$LME_NUMBER) %in% east_atl,]
#get rid of buffer for east atl as well to allow for union
east_atl_spdf_nobuf <- gBuffer(east_atl_spdf, byid=TRUE, width=0)
Spatial object is not projected; GEOS expects planar coordinates
#dissolve polygons into one along coastline
east_atl_spdf_nobuf_agg <- gUnaryUnion(east_atl_spdf_nobuf)Extract bathymetry data from polygon only to make sure we’re limiting to shelf regions above 2000 meters
#crop bathymetry layer to LME subset (continental shelf habitat in LMEs)
raster_extent <- crop(etopo_shelf_raster, extent(east_atl_spdf_nobuf_agg)
#which areas of raster fall within borders?
east_atl_spdf_nobuf_agg_mask <- mask(raster_extent, east_atl_spdf_nobuf_agg)
Error: unexpected symbol in:
"#which areas of raster fall within borders?
east_atl_spdf_nobuf_agg_mask"How to calculate area?
raster::area() Raster objects: Compute the approximate surface area of cells in an unprojected (longitude/latitude) Raster object. It is an approximation because area is computed as the height (latitudinal span) of a cell (which is constant among all cells) times the width (longitudinal span) in the (latitudinal) middle of a cell. The width is smaller at the poleward side than at the equator-ward side of a cell. This variation is greatest near the poles and the values are thus not very precise for very high latitudes. If x is a Raster* object: RasterLayer or RasterBrick. Cell values represent the size of the cell in km2, or the relative size if weights=TRUE
raster::area() SpatialPolygons: Compute the area of the spatial features. Works for both planar and angular (lon/lat) coordinate reference systems. If x is a SpatialPolygons* object: area of each spatial object in squared meters if the CRS is longitude/latitude, or in squared map units (typically meter)
rgeos::gArea Returns the area of the geometry in the units of the current projection. By definition non-[MULTI]POLYGON geometries have an area of 0. The area of a POLYGON is the area of its shell less the area of any holes. Note that this value may be different from the area slot of the Polygons class as this value does not subtract the area of any holes in the geometry.
From wikipedia: earth’s radius is 6,378 km (3,963 mi) at the equator to a near minimum 6,357 km at the poles, with mean of 6378
Now, we will split each coastline raster into latitudinal bins of 1Ëš
Eastern Atlantic
east_atl_spdf_nobuf_mask
#split into northern hemisphere and southern hemisphere
north_extent <- c(xmin(east_atl_spdf_nobuf_mask_1s), xmax(east_atl_spdf_nobuf_mask_1s), 0, ymax(east_atl_spdf_nobuf_mask_1s))
south_extent <- c(xmin(east_atl_spdf_nobuf_mask_1s), xmax(east_atl_spdf_nobuf_mask_1s), ymin(east_atl_spdf_nobuf_mask_1s), 0)
#crop east_atl raster above and below 0
east_atl_spdf_shift_agg_north <- crop(east_atl_spdf_nobuf_mask_1s, extent(north_extent))
east_atl_spdf_shift_agg_south <- crop(east_atl_spdf_nobuf_mask_1s, extent(south_extent))
#all 1Ëš latitude sections for east atlantic
east_atl_north_latitudes <- seq(0, ymax(east_atl_spdf_nobuf_mask_1s), by = 1)
east_atl_south_latitudes <- seq(0, ymin(east_atl_spdf_nobuf_mask_1s), by = -1)
Now, use loop to populate data table with area values
cell_surface_areas <- vector()
for (j in 1:ncell(segment_north)) {
cell <- rasterFromCells(segment_north, j, values = T)
cell_number <- getValues(cell)
cell[] <- segment_north[cell_number]
if(is.na(cell[])) next #skip if value is NA
xmin <- xmin(cell)
xmax <- xmax(cell)
ymin <- ymin(cell)
ymax <- ymax(cell)
area <- grid_cell_area(xmin, xmax, ymin, ymax)
cell_surface_areas <- append(cell_surface_areas, area, after = length(cell_surface_areas))
}